#ifndef CONVECTION_H
#define CONVECTION_H

#include "Embedded/PartialDer.h"

template <int Dim, int pD, int _fitOrder>
class PartialDerOfProd
{
public:
  using iVec = Vec<int, Dim>;
  using iVV  = std::pair<iVec, iVec>;
  using rVec = Vec<Real, Dim>;
  using EmGrid = EmbeddedGrid<Dim>;
  enum { fitOrder = _fitOrder };

  PartialDerOfProd(const EmGrid &aGrid, const Tensor<Real, Dim+1> &aVel);

  void apply(const Tensor<Real, Dim> &aData0, const Tensor<Real, Dim> &aData1, Tensor<Real, Dim> &outData) const;

protected:
  const EmGrid    &emGrid;
  RowSparse<iVec, iVV> matOp;
  const Real      dangerousCN = 1e4; // give a level-1 warning if the condition number is larger than this.
};

template <int Dim, int fitOrder = 5>
class Convection
{
public:
  using iVec = Vec<int, Dim>;
  using rVec = Vec<Real, Dim>;
  using EmGrid = EmbeddedGrid<Dim>;

  Convection(const EmGrid &aGrid, const Tensor<Real, Dim+1> &aVel);

  ~Convection();

  void apply(const Tensor<Real, Dim+1> &aVel, Tensor<Real, Dim+1> &aCvt) const;

protected:
  const EmGrid &emGrid;
  void         *pdpObjs[Dim];
};

#endif //CONVECTION_H
